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Sap transport in trees has long fascinated scientists, and a vast literature exists 
on experimental and modelling studies of trees during the growing season when 
large negative stem pressures are generated by transpiration from leaves. Much 
less attention has been paid to winter months when trees are largely dormant but 
nonetheless continue to exhibit interesting flow behaviour. A prime example is sap 
exudation, which refers to the peculiar ability of sugar maple (Acer saccharum) and 
related species to generate positive stem pressure while in a leafless state. Exper¬ 
iments demonstrate that ambient temperatures must oscillate about the freezing 
point before significantly heightened stem pressures are observed, but the precise 
causes of exudation remain unresolved. The prevailing hypothesis attributes exu¬ 
dation to a physical process combining freeze-thaw and osmosis, which has some 
support from experimental studies but remains a subject of active debate. We ad¬ 
dress this knowledge gap by developing the first mathematical model for exudation, 
while also introducing several essential modifications to this hypothesis. We derive a 
multiscale model consisting of a nonlinear system of differential equations governing 
phase change and transport within wood cells, coupled to a suitably homogenized 
equation for temperature on the macroscale. Numerical simulations yield stem pres¬ 
sures that are consistent with experiments and provide convincing evidence that a 
purely physical mechanism is capable of capturing exudation. 

Keywords: Tree sap exudation; sugar maple; multiphase flow and transport; 
phase change; differential equations; periodic homogenization 


1. Introduction 

The study of tree sap flow has a long history that has given rise over time to the 
concept of the hydraulic architecture of trees [S3]. Despite the extensive literature 
on this subject, several aspects of sap transport remain controversial, including the 
cohesion-tension theory of sap ascent lllSIllSHj; embolism formation and recov¬ 
ery (35] [59], which is ubiquitous in species subject to drought- or freezing-induced 
stresses; and sap exudation in maple and related species such as walnut, butternut 
and birch m- Furthermore, there is a great deal of current interest in the possible 
effects of recent changes in weather patterns on both individual trees and forest 
ecosystems [S] 120], and their connections with sap hydraulics m- The problems 
just described involve complex interactions between sap flow and other phenomena 
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such as nutrient transport, photosynthesis, soil physics, atmospheric dynamics, cell 
growth, etc. Despite the extensive work to date on mathematical and computational 
modelling of trees and their interactions with the environment, many open questions 
remain that can only be addressed by considering sap flow coupled with other pro¬ 
cesses and building integrated models that connect flow and structure at different 
spatial scales and levels of organization 

Sugar maple is a keystone species in the forests of central and eastern North 
America [35] and so is worthy of special attention. Members of the maple family 
are distinguished from other hardwoods by a number of unusual structural and 
functional features that allow them to exude sap during winter [3Il|49j, to gen¬ 
erate unusually high rates of nitrification ED, or to recover from freeze-induced 
embolism US [ST]. The potential impacts of climate change on maple have also at¬ 
tracted recent attention EDIM], motivated by the economic importance of the maple 
syrup industry, not to mention maple’s high timber value. In particular, maple sap 
yields are sensitive to even small variations in temperature or snow cover during the 
harvest season, so that recent unusual weather patterns underscore the importance 
of developing a better understanding of the effects of local environmental conditions 
on sap flow umi]. 

Hundreds of scientific papers have addressed the phenomenon of sap exudation 
during winter when maple trees are leafless and yet still exhibit pressure variations 
that range over 150-180 kPa [31 dn 111119]. However, the precise mechanism driving 
the generation of heightened exudation pressure is still not fully understood [54] . 
The first systematic study appeared in an 1860 article by Sachs [13], who attributed 
exudation pressure to thermal expansion of gas within sapwood or xylem. The next 
major advance in understanding followed from the exhaustive study of Wiegand ES], 
who found to the contrary that thermal expansion of gas, water or wood has min¬ 
imal impact on exudation. Instead, Wiegand proposed a vitalistic or ‘living cell’ 
hypothesis wherein sugar is released into the sap by some cellular activity, which 
gives rise to elevated pressure from osmotic gradients across selectively-permeable 
membranes separating wood cells. Subsequently, this osmotic mechanism figured 
prominently in the literature, although experimental studies have continued to yield 
conflicting results that in turn stimulated development of new theories advocating 
alternate (bio-)physical mechanisms. For example, some authors continued to sup¬ 
port the thermal expansion hypothesis [36] , while others advocated the various roles 
of gas dissolution m. cryostatic suction due to freezing EH], or temperature-induced 
changes in bark thickness [32] . More recent studies have led to a new understanding 
of exudation as a physical process deriving from a combination of freezing and thaw¬ 
ing of sap m with osmosis |50| . Although some experimental evidence supports this 
hypothesis m the precise mechanisms behind the exudation phenomenon is still 
not fully understood. 

We aim to resolve this long-standing open question by developing the first math¬ 
ematical model for the freeze-thaw process in maple. We uncover the essential role 
played by two physical mechanisms whose significance has not yet been recognized 
- namely, root water uptake and freezing point depression due to sap sugar content. 
Using numerical simulations of repeated freeze-thaw cycles, we obtain computed 
exudation pressures that are consistent with experimental results. 

Although the focus of this paper is on developing a complete and physically 
consistent model for sap exudation, our results also have more far-reaching conse- 
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quences. This work affords new insights into the complex multi-physics processes 
occurring in trees and also provides a framework for studying other practical ques¬ 
tions of importance to tree physiologists and maple syrup producers. Our model also 
provides a platform for studying related phenomena such as embolism that occur 
in a much broader range of tree species, as well as evaluating the response of trees 
to changes in environmental variables such as temperature and soil moisture arising 
from various climate change scenarios. 

2. Physical mechanism for sap exudation 

(a) Milburn and OMalley’s hypothesis 

Experimental work up to the 1980’s demonstrated that no single physical mech¬ 
anism is capable of capturing measured winter stem pressures m, and sap exuda¬ 
tion remained an unsolved puzzle until the ground-breaking study of Milburn and 
O’Malley [37]. They proposed a physical mechanism based on freezing and thawing 
of sap, motivated by the unique structural characteristic of xylem in maple (and 
related trees) that a significant proportion of the libriform fibers (or simply fibers) 
are primarily gas-filled rather than being liquid-filled as in most other hardwood 
species [56] . This peculiar feature of the fibers should be contrasted with the two 
other cell types that play an active role in sap transport - vessels and tracheids - 
which are mostly sap-filled and are connected hydraulically to each other via paired 
pits (see figured^). Indeed, recent experiments [TT1|33] suggest that fibers are essen¬ 
tially non-conductive in comparison with other xylem elements because they lack 
end-to-end cell connections and their lateral walls contain mostly unpaired (blind) 
pits that are smaller and fewer in number than the more conductive vessels and 
tracheids. 

Milburn and O’Malley focused on the dynamics of a single fiber-vessel pair as 
pictured in hgure [TJa, and ignored other xylem elements such as parenchyma and 
ray cells. Whereas fibers had previously been thought to play a purely structural 
role, Milburn and O’Malley proposed that as temperature falls below freezing, liquid 
from the vessel is drawn by cryostatic suction through the porous fiber-vessel wall to 
freeze on the inside of the fiber (figured] stages 1-2-3). As a result, any gas contained 
within the fiber is compressed and acts as a pressure reservoir. When temperature 
subsequently rises and ice thaws, the process reverses and the compressed gas bubble 
forces melt-water into the vessel, thereby re-pressurizing the vessel sap (figure |2l 
stages 3-4-1). 

(6) Tyree’s modified hypothesis, with gas dissolution and osmosis 

This freeze-thaw hypothesis was critically evaluated by Tyree m, who proposed 
a modified hypothesis featuring two important additions. First, he recognized that 
gas under pressure will dissolve within an adjacent liquid m and that pressures 
encountered in maple xylem are high enough that any bubbles should dissolve com¬ 
pletely given sufficient time [52]. Therefore, some additional mechanism is required 
to sustain gas bubbles in the fibers. Tyree’s second observation was that measured 
xylem pressures depend strongly on sugar concentration in the vessel sap, 80% of 
which derives from sucrose m, which led him to conclude that sucrose is required 
for exudation. He recognized that although the axial conductivity of fibers is neg- 
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Figure 1: Xylem microstructure, (a) Cross-sectional view of hardwood xylem, 
showing tracheids connected hydraulically to vessels and other tracheids via paired 
pits. Fibers appear similar to tracheids except that they have fewer pits, most of 
which are blind or unpaired. The parenchyma are living cells whose main role is 
carbohydrate storage and so they are ignored here, (b) A fiber-vessel pair ap¬ 
proximated as circular cylinders, showing typical dimensions of the fiber (length 
= 1.0 X 10“^ m and radius = 3.5 x 10“® ) and vessel (L*' = 5.0 x 10“'^ m and 
i?'" = 2.0 X 10“® m). The model domain corresponds to the horizontal cross-section 
through the middle of the diagram. 


ligible in comparison with vessels and tracheids, the lignified cellulose making up 
secondary cell walls should admit a small radial conductivity. He then hypothesized 
that the fiber-vessel wall forms an osmotic barrier that allows water to penetrate 
but not the larger sucrose molecules. Consequently, an additional osmotic pressure 
difference exists between the sweet vessel sap and pure fiber water, which he argued 
is responsible for preventing fiber gas bubbles from completely dissolving. 

Tyree’s modified freeze-thaw hypothesis includes water phase transitions, gas 
dissolution and osmosis, and is currently the prevailing hypothesis for sap exuda¬ 
tion m- It depends strongly on the existence of a hydraulically isolated system of 
fibers and a selectively permeable fiber-vessel wall, both of which have since been 
confirmed experimentally m- Although this evidence is compelling, there has been 
no attempt yet to model this process mathematically (except for a related process 
without phase change in the context of embolism recovery in maple m) and so it 
remains unclear whether this physical description is capable of capturing exudation. 


(c) Three essential physical mechanisms 

Before proceeding further, we extend the freeze-thaw hypothesis just described 
by incorporating three additional mechanisms: 

Gas bubbles in the vessel: Sap (like water) is an incompressible fluid so that in the 
rigid, closed vessel network of a leafless tree there is no mechanism for fiber- 
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Figure 2: Stages in the freeze—thaw cycle. Various stages in the freeze-thaw 
cycle are depicted within an adjacent fiber-vessel pair. Stages 1—>-2^3 depict the 
freezing process: when temperature drops, an ice layer grows on the inner wall of the 
gas-filled fiber as water is drawn via cryostatic suction through the porous cell wall. 
Stages 3—)-4—>-1 depict the reverse process as temperature rises. Note the reversed 
order of phase interfaces inside the fiber between stages 2 and 4. The blue arrows 
denote water transport either through the fiber-vessel wall or between roots and 
vessels. 

vessel mass transfer if vessels are completely saturated with sap. However, the 
existence of gas bubbles within vessels is well-documented in maple [3H1 US] 
and other hardwood species m- Even if xylem pressures were high enough 
to dissolve such bubbles, gas would eventually be forced out of solution upon 
freezing and so at least a transient presence of gas bubbles is unavoidable. 
Therefore, introducing a gas phase in the vessels provides a plausible mecha¬ 
nism for fiber-vessel pressure exchange. 

Sap freezing point depression (FPD): Sap contains dissolved sugars and hence ex¬ 
periences a reduced freezing point compared to pure water according to Blag- 
den’s Law [7], AT/pd = K},Cs/pw, where K], is the cryoscopic constant, Cg 
is sugar concentration, and is water density. For example, sap containing 
3% sucrose by mass experiences a FPD of ATfpd Rs 0.162 °K. Although this 
temperature difference may appear insignificant, we will see that it is actually 
large when considered on the scale of individual cells, and indeed is sufficient to 
account for the existence of ice in fibers while sap in adjacent vessels remains 
in liquid form. This partitioning of ice and liquid in neighbouring fiber-vessel 
pairs induces cryostatic suction that draws liquid out of the vessel to form ice 
on the inner fiber wall. 

Root water uptake during freezing: No previous hypothesis for sap exudation ex¬ 
plicitly considers the role of root water uptake. Furthermore, several studies 
suggest that root pressure in maple has a negligible effect on exudation laisnj- 
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Nonetheless, it is well-known that during winter months trees can draw wa¬ 
ter from the roots if soil temperatures are high enough [JSl [S3], which can 
be caused by an insulating snow cover m- Recent experiments on maple 
saplings (6] [40] have provided the first direct evidence that root water uptake 
occurs in maple during winter while exudation is underway. 

Our aim is now to incorporate these three modifications into a model of the freeze- 
thaw process outlined previously, and then demonstrate that the resulting equations 
are capable of reproducing observed behaviours. 

3. Mathematical formulation 

(a) Outline of the modelling approach 

The freeze-thaw mechanism outlined in the previous section involves processes 
operating on two distinct spatial scales: the microscale corresponding to individ¬ 
ual wood cells with dimensions ranging from 10-100 microns; and the macroscale 
corresponding to the tree stem with diameter tens of centimetres. The derivation 
of our mathematical model for sap exudation therefore divides naturally over these 
two scales. Firstly, we develop microscale equations that capture cell-level processes 
within libriform fibers and vessels, combining the dynamics of freezing, thawing, gas 
dissolution, osmotic pressure, heat transport, and porous flow through the fiber- 
vessel wall. Secondly, we consider heat transport in the entire tree stem and apply 
periodic homogenization to derive an equation for the macroscale temperature that 
incorporates microscale cellular processes via appropriately defined transport coef¬ 
ficients and source terms. We proceed as follows: 

• Start from an existing microscale model for the thawing half of the freeze- 
thaw process [SKIS] in which a 2D periodic microstructure is assembled from 
copies of a reference cell Y containing a single fiber and vessel (see figure 13^). 
The fiber is placed at the centre of Y (where the dashed line denotes the 
fiber-vessel wall) and the remainder of the reference cell corresponds to the 
vessel. This choice of geometry is a mathematical idealization that captures the 
volumes of the fiber and vessel compartments, but is not intended to accurately 
represent the actual layout of wood cells. The governing equations consist of a 
partial differential equation (PDF) for the microscale temperature along with 
five ordinary differential equations (ODEs) for phase interface locations and 
root water volume, coupled nonlinearly through source terms and algebraic 
constitutive relations. 

• Supplement the thawing model with analogous equations for the freezing pro¬ 
cess, which have similar structure but differ slightly depending on the precise 
state of freezing or thawing in the fibers and vessels. 

• Apply periodic homogenization EKig to derive a macroscopic equation for 
temperature that is coupled to the microscale (reference cell) problem at each 
point within the tree stem (see figure [3K)- The macroscopic heat diffusion 
equation contains an integral source term depending on the microscale tem¬ 
perature and capturing all processes on the cellular level. A similar approach 
has been applied to studying protein-mediated transport of water and solutes 
in non-woody plant tissues m- 
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(a) (b) 


Figure 3: Multiscale problem geometry, (a) Idealized microscale fiber-vessel 
geometry, consisting of a square reference cell Y with side length This diagram 
depicts a thawing scenario (stage 4 in figure [2]) wherein a fiber of radius (dashed 
line) contains a gas bubble surrounded by annular layers of ice and liquid water, 
where the gas/ice and ice/water interfaces are concentric circles of radius Sg and s™ 
respectively. The vessel contains a gas bubble (radius r) and liquid sap (water plus 
sugar). The total liquid volume transferred from fiber to vessel is denoted by U. The 
reference cell Y = is divided into two regions separated by a curve T, where 

diffusion on F^ (light blue, outer vessel) is fast and on F^ (dark blue, fiber plus 
fiber-vessel overlap region) is slow, (b) A requirement for homogenization is that 
the tree cross-section can be approximated by a periodic fine-structured domain, 
tiled with copies of the reference cell. The macroscale problem is then solved on a 
homogeneous domain Q having the same size. Radial coordinates on the micro- and 
macroscales are denoted y and x respectively. 

• Exploit radial symmetry on the micro- and macroscales to reduce both PDEs 
for temperature to a single spatial dimension. We will see later in section 3^ 
that the microscale equations need only be solved on the circular sub-region 
F^ in figure [3^ (consisting of the fiber and surrounding vessel overlap region) 
which is clearly radially symmetric. 

(&) Microscale equations for cell-level thawing proeess 

The cell-level model is based on equations already developed for the thawing 
half of the freeze-thaw process by Ceseri and Stockie [8], which were subsequently 
homogenized by Graf and Stockie [H]. We therefore begin by considering an inter¬ 
mediate state in the thawing process corresponding to stage 4 in figure 31 during 
which the vessel sap is completely thawed while the fiber contains both liquid and 
ice. We extend the Ceseri-Stockie model by incorporating additional physical ef¬ 
fects that capture the influence of ice-water surface tension, root water uptake, and 
volume change due to ice/water phase transitions. We discuss some of the most 
important assumptions and modifications next, leaving the reader to consult the 
references [HIIS] for a complete derivation and discussion of assumptions. 

Our model is based on the conceptual diagram in figure [TId that depicts a single 
vessel-fiber pair. Tracheids are not treated separately but instead ‘lumped together’ 
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with vessels because, although they are connected hydraulically to vessels via paired 
pits, they have a much smaller diameter and correspondingly lesser influence on sap 
transport than vessels. Because multiple fibers adjoin and interact hydraulically 
with each vessel, we introduce the parameter N-f representing an average number of 
fibers per vessel, which is estimated from SEM images m as ss 16. Our model 
captures the dynamics of a single fiber and then scales all fiber-vessel flux terms by 
an appropriate factor of N-^. 

We assume that sapwood can be represented as a doubly-periodic array of ide¬ 
alized reference cells Y as pictured in figure [21 where each reference cell contains a 
circular fiber embedded within a surrounding square liquid region representing the 
adjoining vessel. This choice of geometry is made for mathematical convenience in 
the homogenization step, and can be justified because our aim is to derive a sys¬ 
tem of equations that captures the net effect of sap flow and heat transport on the 
microscale, keeping in mind that any specific geometric details will ultimately be 
‘averaged out’ during the homogenization process anyways. 

Our 2D geometry comes with the built-in assumption that axial (vertical) varia¬ 
tions are neglected. In the absence of root water uptake, the model tree behaves as 
a closed system that is essentially in equilibrium. Any pressure differences initiated 
by phase change engender primarily horizontal flow between neighbouring cells, and 
negligible axial flow. Furthermore, we have already shown [8] that phase change on 
the microscale dominates the pressure exchange process and occurs very rapidly (on 
the order of milliseconds). Root water uptake induces an axial flow but this is a 
much slower process; therefore, over the time scales that dominate the microscale 
problem, axial transients may be neglected. 

The fiber is a circular cylinder with length and cross-sectional radius R-^ as 
pictured in figure [3^. Situated at the centre of the fiber is a cylindrical gas bubble 
with time-varying radius Sg(t), outside of which lies an annular ice layer with outer 
radius Siw (t) . The remaining volume extending to the fiber radius contains melt¬ 
water from thawed ice. We note that this configuration is specific to the thawing 
process, and the ordering of ice and water layers would be reversed during freezing. 

The vessel is represented by the portion of the reference cell lying outside the 
fiber-vessel wall (denoted by a dashed line) and the side length I of the reference cell 
is chosen so that the vessel cross-sectional area equals that of a cylinder of radius 
R". Keeping in mind that there are actually fibers connected to each vessel, we 
require that (. satisfy the area constraint 


£2 = 7r7V'^(i?/)2 +7r(i?«)2. 


(3.1) 


Within the vessel is a gas bubble of radius r{t), which is surrounded by liquid sap 
owing to the FPD effect that lowers freezing temperature below that in the fiber. 
The cumulative volume of melt-water flowing through the porous fiber-vessel wall 
is denoted by U{t) and is measured positive from fiber to vessel. The final variable 
that determines the local state of the fiber-vessel system is the volume of root water 
uptake, denoted Urootft). 

We may now formulate a first-order system of five ODEs describing the time 
evolution of Siw, Sg, r, U and Uroot- The fiber ice-water interface is governed by the 
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Stefan condition [DIS] 




{Ew ~ Ei) 


VT 2 • n + 


dtU 

27r Si'uj Ld 


(3.2) 


where VT 2 ■ ft represents the normal temperature derivative on the interface (i.e., 
the curve along which temperature equals the melting (or freezing) point T^) and 
the final term accounts for the volume of water transferred between fiber and vessel. 
This form of the Stefan condition assumes that liquid motion induced by phase 
density differences is negligible [1]. The microscale temperature T 2 {y,t) is obtained 
as the solution of a heat diffusion equation that will be stated in the next section, 
where the microscale spatial coordinate is y. The parameters and kw denote 
density and thermal conductivity of liquid water, while {E^j — Ei) is the enthalpy 
difference between water and ice (also called the latent heat or enthalpy of fusion) 
at locations where T 2 = T^. The effects of thermal expansion are known to be 
relatively small [56] and so have been neglected here. 

Imposing mass conservation yields an equation for the fiber gas bubble radius 
(which in this thawing scenario is a gas-ice interface) 


{Pw Pi)^iwflt^iw PwfdtU 

where pi is the density of ice. An equation for the vessel gas bubble radius follows 
from a similar mass conservation argument 



dtr = 


N^dtU + dtUroot 
2TTrL'’ 


(3.4) 


where L" denotes the length of a vessel. This last equation expresses the balance 
between water flux from neighbouring fibers and the slight volume change stemming 
from the water/ice density difference. The effect of gas dissolution has been omitted 
here but will be incorporated below in the gas density; this approximation was al¬ 
ready justified in (2, which showed that incorporating dissolution in these equations 
has negligible impact on the bubble radii Sg and r. 

Darcy’s law governs liquid water flux through the porous fiber-vessel wall 


dtU = - 


Nf L' 


Posm +p{{t) 


(3.5) 


where the wall is characterized by hydraulic conductivity Af’ and surface area A. The 
pressure term in square parentheses derives from four contributions: liquid pressure 
in the vessel {p"^) and fiber (p{,), osmotic pressure (posm), and capillary pressure 
{p{) due to ice-water surface tension [15]. This latter contribution, also known as 
cryostatic suction, follows hand-in-hand with FPD and arises whenever ice lies on 
the inside surface of the wall and liquid sap is present on the vessel side, since then 
the small capillary pores in the adjoining wall (with radius reap) contain both ice 
and liquid. For the thawing scenario under consideration here, water lies on both 
sides of the fiber-vessel wall and so p{ = 0; however, other stages in the freeze-thaw 
process can give rise to non-zero p{ as detailed in the next section (see also figure]?]). 

The final ODE comes from another application of Darcy’s law to root flux 

dtUroot = max I - (p((,(t) -psoti), o|, (3.6) 
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where is the root hydraulic conductivity and Ar denotes the portion of root 
surface area corresponding to a single vessel. The cut-off function ‘max{ •, 0}’ ensures 
that water only flows inward from soil to roots and not outward, which is consistent 
with experiments that demonstrate root outflow can be a factor of five smaller than 
that for inflow [21] . Indeed, studies of root water transport in a variety of tree species 
show that root conductivity can vary with factors such as temperature m , root 
age m], and time of day m or season [33]. Many authors attribute this selective 
control of water transport to membrane proteins known as aquaporins m- 

In the preceding discussion we introduced a number of constant parameters 
whose values are listed in table |TJ The remaining symbols correspond to interme¬ 
diate variables whose definitions we provide next. First, the density of gas in the 
fiber and vessel bubbles depends on initial values of density and volume, modified 
to account for dissolved gas according to 



/ F/(0)-bgFi(0) \ 

V V/ + HVd ) 




Pa 


f V^{Q) + HVZ{Q) \ 
V V- + HV- ) 


pm. 


(3.7) 


where H is the dimensionless Henry’s constant for air in water. The various 
volumes are determined from the cylindrical cell geometry as 


■vJ2 


V/=TTLfsl, Vg^=TTrr 

V-f = nLf ((i?0' - sL) , v: = ttL^ ((E’’)" - . 


phase 


(3.8) 

(3.9) 


The corresponding gas pressures are given by the ideal gas law as 

P^^^ _ Pa^^ 


pi = 


Mn 


Pa = 


Mr, 


(3.10) 


where is the universal gas constant and Mg is the molar mass of air. The water 
and gas pressures in both fiber and vessel differ by an amount equal to the capillary 
pressure, which is determined by the Young-Laplace equation as 


f f 2(t„u, 

Pin=Pi -: 

^ c 


2a. 


Pw = Pa 


gw 


(3.11) 


where ag^, is the air-water surface tension. The osmotic pressure across the fiber- 
vessel wall depends on sap sugar concentration according to 


Posni = Cs^T. 


(3.12) 


Finally, the sap sugar content induces a reduction in freezing temperature that obeys 


Tm,sap =Tm — ^Tfpd = — 


KbCs 


(3.13) 


(c) Equations for other phase transitions 

In the previous section we developed equations specific to the thawing process, 
during which the vessel is completely thawed and the fiber contains a mix of gas, 
water and ice (see stage 4 in figure |3|). We describe next how these equations should 
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Table 1: Model parameters for base case simulation 

(Unless cited otherwise, all parameter values are taken from [S]) 


Symbol 

Description 

Values 

Units 

Microscale variables (functions of time t and space x^y): 
Siwi Sg interface locations in fiber 


m 

r 

vessel bubble radius 


m 

u 

water volume flowing from fiber to vessel 


m® 

Uroot 

root water volume uptake 


m® 

V 

volume 


m® 

p 

pressure 


Pa 

p 

density 


kg m~^ 

Subscripts: i, w, g for ice, water/sap, gas 

Superscripts: /, v for fiber, vessel 

Tree structural parameters: 

A area of fiber—vessel wall 

6.28 X 10“® 


Af 

root area for a single vessel |15| 

1.14 X 10"® 


i 

side length of reference cell, equation (13.111 

4.33 X 10“® 

m 

Lf 

length of fiber 

1.0 X 10“® 

m 

L” 

length of vessel element 

5.0 X 10-"^ 

m 


conductivity of fiber-vessel wall 

5.54 X 10-1® 

m s“^ Pa~^ 


conductivity of roots |47l I53| 

2.7 X 10-1® 

m s“^ Pa~^ 

AT/ 

number of fibers per vessel 

16 

- 

Rt 

inside radius of fiber 

3.5 X 10-® 

m 

K“ 

inside radius of vessel 

2.0 X 10-® 

m 

't' cap 

radius of pores in fiber-vessel wall |28| 

2.80 X 10-1 

m 

w 

thickness of fiber-vessel wall 

3.64 X 10-® 

m 

Water phas 

e properties: 

ice, liquid 


Ci , C-uj 

specific heat capacity 

2100, 4180 

J °K-i kg-i 

Eii Ew 

enthalpy at Tm 

574, 907 

kj kg-i 

1 l^w 

thermal conductivity 

2.22, 0.556 

W m-i °K-i 

pi ) Pw 

density 

917, 1000 

kg m~^ 

^ iw 1 ^ gw 

surface tension |18| 

0.033, 0.076 

N m-l 

Coo 

regularization parameter, equation II3.32D 

1.0 X 10i 

J °K-l kg-l 

Physical constants: 

H Henry’s constant for air in water 

0.0274 


Kb 

cryoscopic (Blagden) constant 

1.853 

kg °K mol-l 

Mg 

molar mass of gas (air) 

0.029 

kg mol“^ 


universal gas constant 

8.314 

J °K-l mol-l 

T-ni. 

melting point for pure water 

273.150 

°K 

‘Base case’ 

Cs 

simulation: 

sap sugar concentration (3% by mass) 

87.6 

mol m~^ 

Psoil 

soil pressure at roots = p^(^) 

2.03 X 10® 

Pa 

R 

tree cross-sectional radius 

0.035 

m 

Ta{t) 

ambient temperature 

[-10,20] +T^ 

°K 

Tm,sap 

melting point for sap = Tm — Kf^Csjpw 

272.988 

°K 


be modified to capture other freeze-thaw states in the fiber and vessel. In particular, 
we account for the fact that phase interfaces can appear or disappear whenever ice 
completely thaws (or liquid completely freezes), as well as the reversal of the ice and 
water layers in the fiber during freezing and thawing. 

In fact, many equations remain unchanged throughout the entire freeze-thaw 
cycle, with the exception being those for Sg , s™ and pf. The required modifications 
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Figure 4: Microscale equations for all stages of the freeze thaw process. 


for each case are listed in figure 01 referenced by the numbered stages in figured] We 
emphasize that the ice-water interface lies within pores in the fiber-vessel wall and 
forms a mushy layer wherein both solid and liquid phases coexist in the pore space. 
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This type of phase interface (called a frozen fringe in the context of ice lensing in 
soils m) differs from an idealized gas-water interface in that the interfacial pressure 
jump p{ increases with ice volume fraction according to 

f _ 2ui^ V/ 

“ reap V/ + vi ’ 

where ai^ represents the ice-water surface tension and are the corresponding 
volume fractions. Finally, we note that when both fiber and vessel are completely 
frozen (stage 3) the equations for r, U and Uroot also drop out of the system. 

(d) Homogenized equation for temperature 

The equations derived in the preceding two sections govern microscale processes 
at the cellular level whereas on the macroscale the temperature is of primary interest, 
and it is transport of heat between the external (ambient) environment and the 
interior of the tree stem that drives the freeze-thaw process. Clearly, there exists a 
two-way interaction between the global temperature and the local fiber-vessel state, 
wherein temperature governs phase change dynamics in fibers and vessels, while 
cellular processes in turn influence heat transport through the Stefan condition and 
local phase volume fractions. To simplify this complex multiscale problem, we exploit 
a separation in spatial scales reflected in the fact that state variables describing the 
fiber-vessel configuration are essentially ‘invisible’ on the macroscale except through 
their effect on heat transport properties of the sap- and gas-filled wood. 

Because of the repeating microstructure of wood, this problem is ideally suited 
to the application of periodic homogenization. The philosophy behind this approach 
is to solve at each point in space a local problem on a reference cell Y that deter¬ 
mines the solution state on the microscale. By using an appropriate homogenization 
or averaging procedure, the effect of microscale variables on the macroscale may 
then be incorporated into equations for the global solution variables. One technical 
requirement is that the reference cell must divide into two sub-regions, Y = UY^, 
separated according to whether heat diffusion is fast (in Y^, the outer portion of 
the vessel) or slow (in an overlap region covering the fiber and the remainder 
of the vessel). The result is two heat equations: one governing temperature on the 
macroscopic domain O and the second on Y^ x O. When these two equations are 
coupled together, we obtain a two-scale temperature solution on the domain Y x fl. 
Instead of fully coupling the micro- and macroscale equations, this homogenization 
approach leads naturally to a simpler system of equations that captures the essential 
aspects of coupling between scales. A similar homogenization approach has been ap¬ 
plied by Chavarria-Krauser and Ptashnyk to a model of water and solute transport 
in plants m- 

The dynamics of heat transport are best described using a mixed formulation 
written in terms of both temperature and specific enthalpy, which are denoted re¬ 
spectively by T 2 {x,y,t) and E 2 {x,y,t) on the reference cell region and Ti{x,t) 
and Ei(x,t) on the macroscale. The variables Ti and Ei depend on time t and the 
global spatial coordinate x, whereas microscale quantities have an additional de¬ 
pendence on the reference cell Y through a local spatial coordinate y. Temperature 
and enthalpy are not independent variables but instead are related via the piecewise 
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linear function 

r j-E, iiE<E,-6,, 

T{E) = < Tm + Ei — 5i < E < Eyj + 5yj, (3.32) 

[_ + '^{E — E,„), if ifu, + 5.U) ^ E. 

We introduce the large parameter Coo (taking Coo = 10^ in practice) to impose a 
small but nonzero slope (1/coo) on the central plateau region where T « Tm- We 
also make use of the fact that Ei = CiTm and choose 


2(Coo Cj) 


and Sw 





(3.33) 


so that the function T{E) is continuous. This form of T{E) is a regularization of the 
exact temperature-enthalpy relationship |55j that avoids numerical instabilities in 
the calculation of temperature and also recovers the exact (piecewise linear) result 
in the limit as Coo —t oo and <5^, —>■ 0. 

During the homogenization procedure |19| . we find that heat transport in the 
reference cell must only be treated on the sub-region where temperature obeys 


Cn,dtT 2 - Vy ■ {D{E 2 )VyT 2 ) = 0 in t) x D, (3.34) 


and D{E 2 ) is a thermal diffusion coefficient that is a piecewise linear and continuous 
function of enthalpy m 

HE < Ei, 

l), HE,<E<E,„ (3.35) 

if E,t, < E. 

We employ this nonstandard definition of D (instead of the usual thermal diffusivity 
having units m^ s“^) so that we can factor out the specific heat, thereby allowing 
the same coefficient to be used in both this microscale heat equation and the mixed 
temperature-enthalpy form we develop below for the macroscale. We include an 
explicit time- and global space-dependence in Y^{x,t) to emphasize the fact that 
the ice region within the fiber is bounded by a moving water-ice interface, and that 
the fiber configuration varies from point to point throughout the tree stem. On the 
water-ice interface (corresponding to the inner boundary of T^) the temperature 
equals the melting point value 

T 2 = Tm ondY^ix,t) xn. (3.36) 


D{E)={ ^ 


Pi 

E-Ej 

E.^-Ei 

k-a 

Pu 


\ pw 


We thereby obtain the macroscale temperature equation 

dtEi - V, • (nD(£;i)V,Ti) = ^ f D{E2)VyT2 ■ n dS in D, (3.37) 

D \ Jr 

where the coupling with microscale variables is embodied in a surface integral term. 
The factor If multiplying the diffusion coefficient is a 2 x 2 matrix whose entries 
depend on the reference cell geometry according to 

^ij ~ |yi I J ^ E ’S/y^i) Ay, (3.38) 
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for i,j = 1,2. Here, 5ij is the Kronecker delta symbol and are solutions of a 
standard reference cell problem on [5]. The temperature on the outer surface of 
the tree is held at the ambient value 

Ti = Ta{t) on dn. (3.39) 

Finally, the micro- and macroscale solutions are coupled by matching temperature 
on the interior boundary 

T 2 = Ti on r X fl. (3.40) 

In summary, the governing equations consist of a system of differential-algebraic 
equations (I3.2I) - (I3.13I) and (I3.34I) - (I3.36I) for the microscale temperature and fiber- 
vessel state variables within each local reference cell. These are supplemented by 
equations (I3.37I1 - (I3.40I) for the macroscale temperature on H. Both problems are 
solved at each spatial point x G and the two solutions are coupled by means of 
the integral source term in (I3.37P and the boundary condition p.40l) . The geometry 
of the local reference cell is also incorporated into the macroscale problem via the 
(constant) pre-factors H multiplying the diffusion coefficient in (13.371) . 

4. Simulating daily freeze—thaw cycles 

(a) Numerical solution algorithm 

The radial symmetry of both micro- and macroscale domains implies that all 
solution variables can be written as functions of a single radial coordinate and time. 
We use a method of lines approach and discretize the temperature variables in 
space using finite elements, yielding a large system of time-dependent ODEs. When 
combined with the ODEs and algebraic equations governing microscale fiber-vessel 
dynamics, the resulting coupled system is integrated in time using a standard ODE 
solver. The spatial discretization on the two scales proceeds as follows: 

• Microscale (cell-level) equations: The fiber ice temperature is assumed to be 
a uniform 0 °C, and gas temperature is also taken constant since the thermal 
diffusivity of gas is so much larger than that for either ice or water. Therefore, 
the PDE (13.341) for temperature on must only be solved on the annular 
region between T and the phase interface s™ (see figure [3^). We find that 
sufficient accuracy is obtained for T 2 by using only 4 radial grid points within 
the annulus. Because the phase interface evolves in time, we use a moving 
mesh approach wherein the motion of grid points introduces an additional 
‘grid advection’ term that is proportional to the mesh point velocity [33] . 

• Macroscale (tree-level) equation: The tree stem is similarly divided into equally- 
spaced radial points, and here we find that taking 20 grid points yields suffi¬ 
cient accuracy in Ti. Owing to radial symmetry, the integral source term in 
(13.371) reduces to multiplication by the curve length |r|. The factors H de¬ 
pend only on the reference cell geometry and so can be pre-computed at the 
beginning of a simulation. 

We employ an efficient split-step approach where in each time step the reference cell 
problem is solved for the microscale variables, and then the macroscale temperature 
equation is solved by holding the microscale variables constant. 
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The algorithm described above has been implemented in Matlab using the built- 
in stiff ODE solver odel5s to integrate the equations in time. The only algorithmic 
detail remaining to be described is the switching between equations required as 
phase interfaces appear or disappear. We can capture this switching simply and 
robustly using the Events option provided in the ODE solver suite, which signals 
an event based on zero-crossings of an ‘indicator function’. During any portion of 
the freeze-thaw cycle, the indicator function is set equal to either the thickness of 
a phase interface or the difference between the phase temperature and the melting 
temperature. When the indicator crosses zero, the time integration halts, equations 
are modified appropriately, and the ODE solver is restarted using the new set of 
equations and taking the current solution as the new initial state. The time in¬ 
tegration then proceeds until the next phase change event is signalled. A typical 
simulation covering 4 daily temperature cycles requires between 30-45 minutes of 
clock time on an Apple MacBook Pro with 2.3 GHz quad-core Intel i7 processor. 

(b) Choice of parameters 

The algorithm just described is used to simulate freeze-thaw dynamics in a 
typical base case scenario for which all parameters are listed in table [TJ We take 
a ‘sapling’ of diameter 0.07 m consisting entirely of sapwood. The sugar content 
of maple sap ranges from 1-5% by mass [49] and so we choose a representative 
value of 3% that induces a vessel FPD of = 0.162 °C. To mimic temperature 

variations during late winter, we let ambient temperature vary sinusoidally between 
— 10 and -1-20 over a 24-hour period (this range is somewhat extreme but is 
chosen to correspond with the experiments of Ameglio et al. [^ that we will describe 
shortly). We begin with a freezing event and initialize the tree in a thawed state 
with uniform temperature 0.35 “’C, just slightly above the freezing point. Each fiber 
initially contains gas and water with 75% gas by volume, whereas the vessel has a 
much smaller initial gas content of 8%. 

There remain two parameters whose values we have not been able to obtain rea¬ 
sonable estimates from the literature - root hydraulic conductivity and capillary 
pore radius reap - and so we have had to adjust their values in order to match numer¬ 
ical results with experimental data. First, we choose = 2.7 x 10“^® m s“^ Pa“^ 
so that pressure and root uptake vary over time scales similar to those observed 
in experiments min]. Then, we take reap = 2.8 x 10 ^ m so that the exudation 
pressure build-up is within the observed range of 80 to 150 kPa [nma. This pore 
size is also consistent with that measured in other membranes that hinder transport 
of sucrose molecules [55]. 

(c) Base case: Pressure build-up during temperature eycling 

Using these base case parameters and initial conditions, we perform two nu¬ 
merical simulations: one with root water uptake corresponding to a soil pressure of 
Psoii = 203 kPa, and a second with no root uptake (e.g., consistent with a completely 
frozen soil). Vessel sap pressures are compared in figure [5^, and in both cases we 
observe a periodic variation in pressure synchronized with daily temperature fluctu¬ 
ations. Without root uptake, the vessel pressure simply oscillates between two fixed 
values of 20 and 200 kPa and there is no pressure build-up over multiple freeze-thaw 
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Figure 5: Comparison of base case simulation with Ameglio’s experi¬ 
ments. (a) Simulated vessel sap pressure (middle, with and without root water) 
and cumulative root water uptake (bottom) in response to an imposed periodic am¬ 
bient temperature (top). The vertical dotted lines highlight times when temperature 
crosses the freezing point. Two primary features of the vessel pressure curve are the 
amplitude of pressure oscillations in each daily cycle (Api, arising from ice-water 
capillary effects) and the residual pressure increase at the end of a cycle (Ap 2 , due 
to root water uptake), (b) Ameglio et aids experiments on black walnut [3] (repro¬ 
duced with permission of Oxford University Press). The measurements relevant to 
our study are ‘P control’ (sap pressure) and ‘T trunk’ (temperature). 


cycles. However, when root uptake is included there is a gradual pressure increase 
superimposed on the background oscillations, with a total increase (measured from 
the local maximum in each cycle) of roughly 80 kPa over the four days. The accom¬ 
panying plot of total root uptake in figure!^ shows that the majority of root water 
is absorbed during the first freeze-thaw cycle, followed by a more gradual uptake 
that is essentially complete after 3 days. 

We next draw a direct comparison with the experiments of Ameglio et al. [5] 
who studied black walnut trees {Juglans nigra) in a controlled laboratory setting 
where the living stump of an excised tree branch was connected via a sealed pipe 
to a pressure transducer. We calculate vessel sap pressure in our simulations as an 
average pressure across the stem cross-section to be as close as possible to such a 
transducer measurement. We are unaware of any comparable data for sugar maple, 
but we claim that a meaningful comparison may still be drawn with Ameglio’s 
results since black walnut is closely related to maple and undergoes exudation under 
similar conditions [21 HU US]. The curves to focus on in Ameglio’s figure [5Jd are the 
air temperature (labelled ‘T trunk’) and stem pressure (labelled ‘P control’). 

The qualitative agreement between simulated and experimental pressures is re¬ 
markable considering the complexity of the processes involved and the minimal 
parameter fitting required. The overall shape of pressure curves is similar, with each 
freeze-thaw cycle exhibiting a rapid increase whenever temperature exceeds the 
freezing point. The pressure then attains a maximum, after which there is a slight 
decrease over roughly 6-8 hours, followed by a rapid drop as ambient temperature 
crosses the freezing point again. We remark that there is also a rough quantita¬ 
tive agreement between simulations and experiments in that pressure oscillations 
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have an amplitude of 80 to 100 kPa, and the total pressure build-up over four days 
also is similar. On the other hand, the maximum value of our simulated pressure is 
290 kPa, which is almost double the 160 kPa observed in the black walnut experi¬ 
ments; however, it is possible that more time is needed for the experiment to reach 
steady state, not to mention that there are species-specific differences that could 
influence pressure. 

A more quantitative comparison can be drawn based on two characteristic fea¬ 
tures of the pressure in figure [5^ labelled as Api and Ap 2 ■ The first corresponds to 
the amplitude of oscillations in the absence of root uptake, which derives mainly from 
cryostatic suction and so can be estimated using the formula Api 2aiwfrcap ~ 
236 kPa. This value is close to the computed amplitude of the ‘no root uptake’ curve, 
as well as to the rise in vessel pressure during the initial thawing event for the ‘base 
case’. The second feature Ap 2 captures the exudation pressure build-up during the 
first freeze-thaw cycle which arises mainly from root water uptake. Because this 
additional water acts to compress the gas in fiber and vessel, we apply the differ¬ 
ential form of the ideal gas law at constant temperature, Ap 2 ~ —pAVfV, during 
the first freezing event. Substituting values of p ~ 200 kPa for the initial vessel 
pressure, AV 0.4 cm^ for the root water volume uptake (taken from figure [S^) 
and V ~ 1.15 cm^ for the initial gas volume in a slice through the tree cross-section 
(with thickness equal to that of the reference cell, L-f), we obtain |Ap 2 | ~ 69 kPa. 
The correspondence between this estimate and the computed value of 50 kPa is 
reasonable, considering that it ignores effects such as gas dissolution. 

Despite the abundance of experimental data available for sugar maple nil HI US 
iS], most experiments measure sap outflux from tapped trees m rather than the 
‘closed system’ corresponding to an untapped tree that we consider here. Other mea¬ 
surements have been taken of excised wood samples rather than living trees, while 
yet others were taken in uncontrolled external conditions with irregular variations in 
pressure and ambient temperature. Consequently, we hesitate to attempt a detailed 
comparison between any of these experiments and our simulations; nonetheless, we 
can still draw a few quantitative comparisons. For instance, Tyree m performed 
experiments on excised maple branches that absorbed water at a maximum rate 
of 12 cm^/h; for similar sized branches, our model yields a comparable maximum 
absorption rate of roughly 13 cm^/h as well as qualitatively similar solution pro¬ 
files. Another experiment by Johnson et al. [3^ yielded total root uptake of 2.0 cm^ 
during freezing, followed by a much smaller uptake of 0.1 cm^ during a subsequent 
freezing event. We see similar qualitative behaviour in our simulations, as well as 
measuring 2.2 and 0.2 cm^ of water absorbed during the first and second freeze, 
respectively. 


(d) Two crucial mechanisms: Root water uptake and FPD 

To evaluate the relative importance of the various physical mechanisms, we 
present in figure |6] the base case pressure and root water uptake alongside simula¬ 
tions with each of the following mechanisms ‘turned off’: FPD, root uptake (repeated 
from figure [S^) , osmosis and gas dissolution. The first two mechanisms clearly have 
the greatest impact on the build-up of exudation pressure. We already discussed 
the crucial role of root uptake in facilitating pressure accumulation over multiple 
freeze-thaw cycles. This effect is underscored by the plots in figure [7^ depicting the 
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(a) 


(b) 



time [d] 



Figure 6: Comparison of various physical mechanisms, (a) An investigation 
of the relative importance of various physical effects, depicting pressure with each 
of the following mechanisms turned off: FPD, root water uptake, osmosis, gas dis¬ 
solution. The ‘base case’ and ‘no root uptake’ curves are repeated from figure [S^ for 
easy comparison, (b) Corresponding curves for root uptake. 


(a) 



time [d] 


(b) 



time [d] 


Figure 7: Sensitivity of exudation pressure to parameters, (a) Root hy¬ 
draulic conductivity, in m s“^ Pa“^. (b) Sugar content in %. 


pressure response when root conductivity varies between zero (no uptake) and 
nearly ten times the base value. 

Without FPD, the vessel pressure remains nearly constant and there is minimal 
root uptake, whereas without osmosis the vessel pressure increases. We therefore 
conclude that the predominant impact of sugar on exudation is through FPD rather 
than osmosis, and even though ATfpd is small it nonetheless plays a critical role in 
facilitating pressure transfer between fiber and vessel. This dependence is illustrated 
further by figure [71d, where sugar content is varied between 0 and 7% and we ob¬ 
serve that both net pressure build-up and oscillation amplitude increase with sugar 
content. 

One assumption requiring further investigation is that of zero conductivity to 
root outflow in (EH), which we motivated by citing experimental results that ex¬ 
hibit a small but still nonzero root outflow m- To study this outflow effect, we 
take four different outflow conductivities equal to the inflow value scaled by a 
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(a) 


(b) 



time [d] 



Figure 8: Sensitivity to root outflow, (a) Pressure curves with non-zero con¬ 
ductivity to root outflow, where the ratio of outflow-to-inflow ratio varies between 
0 and 1. (b) Corresponding plots of root uptake. 


factor between 0 and 1 (where 0 corresponds to the base case). The results in fig¬ 
ure [8] clearly show that allowing even a small outflow has a major influence on the 
root water uptake by preventing accumulation of water over multiple freeze-thaw 
cycles and thereby reducing build-up of exudation pressure. Because of the obvious 
sensitivity of these results to root outflow, a more extensive experimental study of 
root conductivity in maple is warranted. 

We end this section by addressing the seemingly counter-intuitive result in fig- 
ureE^ that introducing osmosis decreases vessel sap pressure. This result can be most 
easily explained by considering the water flux equation (13.51) over a long enough time 
that the fiber and vessel have reached a quasi-steady state and dtU « 0. Then (13.51) 
reduces to the simple pressure balance 




- Posm + p: 


f 


0 . 


The ice-water capillary pressure p{ is a constant, and our simulations show that 
osmosis has relatively small impact on fiber bubble size and pressure (the latter 
effect was discussed in 0). Therefore, the primary influence of osmosis is within 
the vessel: osmotically-driven flow from fiber to vessel compresses the vessel bubble 
which not only increases the vessel gas pressure Pg, but also increases the capillary 
pressure term (via a reduction in bubble radius r). The contribution from surface 
tension dominates and so the net effect is actually a decrease in vessel sap pressure 
p^, which is consistent with figure and the results reported in [5]. 


(e) Phase change dynamics on the microscale 

When a completely frozen tree warms above 0 °C during the day, a thawing 
front develops near the bark (wherein water and sap are frozen ahead of the front 
and thawed behind) and advances into the stem; an analogous scenario occurs upon 
freezing. Clearly, the ‘interesting’ solution dynamics will occur in the vicinity of 
this front, and hence knowledge of phase change on the microscale is desirable for 
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(a) (b) 




0.005 0.01 0.015 0.02 0.025 


(c) (d) 




Figure 9: Local phase change dynamics. Plots of temperature versus radius 
and enthalpy for a fixed time in the middle of a freezing event (a,b, top) and a 
thawing event (c,d, bottom). Points correspond to discrete solution values on an 
equally-spaced radial grid and are coloured according to the current state of hber 
and vessel: blue if both are frozen; red if both are thawed; purple if hber is frozen 
and vessel is thawed. 


understanding solution behaviour. Over a century ago, Wiegand |56j recognized the 
existence of freezing and thawing fronts that ‘penetrate the wood in a wave-like 
manner’ and in which ‘but few cells would actually take part in the production of 
pressure at any one time’; however, there has so far been no attempt to develop a 
mathematical model for this phenomenon. In particular, the role of FPD in govern¬ 
ing the progress of these phase transitions throughout the sapwood has not been 
investigated before. 

Phase change dynamics are most easily studied by means of a temperature- 
enthalpy diagram as depicted in figures which are each taken at a fixed time 
during a freeze or thaw event. Both plots feature a plateau region at the melting 
temperature, which has a horizontal extent equal to the enthalpy of fusion. Note 
that there are two distinct melting temperatures in fiber and vessel equal to and 
Tm,sap = Tm — ^Tfpd respectively. The corresponding plots of temperature versus 
radius are shown in figures IHbjd which depict the local state of each point within the 
tree stem. For example, in the freezing case (top) the three grid points closest to the 
stem centre are completely thawed (red), the outermost point is frozen (blue), and 
the intervening points are undergoing freezing (purple). Owing to FPD, water in the 
fiber freezes before the vessel sap, thereby introducing a time delay in formation of 
ice between the fiber and vessel. 

For both freezing and thawing, the bulk of the stem is in a state located at the 
leading edge of the enthalpy plateau (right edge for freezing, left edge for thawing). 
This behaviour can be explained by considering the local rate of phase change: 
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conservation of energy at a phase interface is expressed mathematically using the 
well-known Stefan condition, which states that the rate of freezing (or thawing) is 
proportional to the temperature gradient. Referring to figures[5)3,d, the temperature 
difference between adjacent points is smaller near the tree centre and larger near the 
bark. Consequently, at any location in the tree a freezing event begins within the 
fiber as a slow process, followed at a later time in the vessel which freezes relatively 
quickly. In contrast, a thawing event begins with a slow thawing of the vessel sap, 
followed by rapid thawing in the fiber. 


5. Concluding remarks 

We have developed the first complete mathematical model for the tree sap exuda¬ 
tion process based on a prevailing freeze-thaw hypothesis. We introduced a number 
of additions to this hypothesis, and identified root water uptake and freezing point 
depression (FPD) as the two main driving mechanisms for sap exudation. In par¬ 
ticular, we showed that the primary mechanism whereby sugar induces exudation 
pressure is via the FPD and not osmosis as was previously believed. Numerical sim¬ 
ulations of the governing equations demonstrate qualitative and quantitative agree¬ 
ment with experimental data on sugar maple and the related species black walnut. 
The quality of agreement is striking considering that the model parameters were 
determined using a minimum of parameter fitting. Our work clearly demonstrates 
the need for further experiments on sugar maple that parallel the work of Ameglio 
et al. on walnut [3]. Our model results lead to the important conclusion that FPD 
is a primary driver of sap exudation, which also requires experimental validation. 
Furthermore, because we have only rough estimates at present for two of the model 
inputs - capillary pore size reap and root conductivity J/fr (especially differences 
between conductivity to inflow and outflow) - more accurate measurements of these 
parameters are also required. 

Our model provides an ideal platform from which to investigate other problems 
related to sap flow in maple and related species. First of all, we aim to extend 
our current model of a 2D stem cross-section to three dimensions. This will permit 
us to incorporate variations in gravitational pressure head and sugar concentration 
with height [56] and to study problems of practical importance to the maple syrup 
industry such as optimizing tap-hole placement or determining sensitivity to changes 
in soil or climatic conditions. Finally, there are a number of intriguing parallels 
between exudation and the phenomenon of freeze-induced winter embolism 113117] 
that are also worthy of future investigation. 
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